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Abstract 


A method is proposed for determining the motion of a body relative to a fixed 
environment using the changing image seen by a camera attached to the body. The 
optical flow in the image plane is the input, -while the instantaneous rotation and trans¬ 
lation of the body are the output. If optical flow could be determined precisely, it would 
only have to be known at a few places to compute the parameters of the motion. In 
practice, however, the measured optical flow will be somewhat inaccurate. It is there¬ 
fore advantageous to consider methods which use as much of the available information 
as possible. Wc employ a least-squares approach which minimizes some measure of the 
discrepancy between the measured flow and that predicted from the computed motion 
parameters. Several different error norms are investigated. In general, our algorithm 
leads to a system of nonlinear equations from which the motion parameters may be 
computed numerically. However, in the special cases where the motion of the camera is 
purely translational or purely rotational, use of the appropriate norm leads to a system 
of equations from which these parameters can be determined in closed form. 


This report describes research done at the Artificial Intelligence Laboratory of 
the Massachusetts Institute of Technology. Support for the laboratory’s artificial 
intelligence research is provided in part by the Advanced Research Projects Agency of 
the Department of Defense under Office of Naval Research contract N00014-80-C-0505 
and i.u part by National Science Foundation Grant MCS77-07569. 

© Massachusetts institute of technokkx kh 











1. Introduction 

In this paper we investigate the problem of passive navigation using optical flow 
information. Suppose we are viewing a film. We wish to determine the motion of the 
camera from the sequence of images, assuming that the instantaneous velocity of the 
brightness patterns, also called the optical How , is known at each point in the image. 
Several schemata for computing optical flow have been suggested (e.g. [2], 3 , 5]). 
Other papers (e.g. [9], [10], [11]) have previously addressed the problem of passive 
navigation. Three approaches can be taken towards a solution which we term the 
discrete, the differential and the continuous approach. 

In the discrete approach, information about the movement of brightness patterns 
at only a few points is used to determine the motion of the camera. In particular, using 
such an approach, one attempts to identify and match discrete points in a sequence of 
images. Of interest in this case is the photogrammetric problem of determining what 
the minimum number of points is from which the motion can be calculated for a given 
number of images [10], [11], [14]. This approach requires that one tracks features, or 
identifies corresponding features in images taken at different times. 

In the differential approach, the first and second spatial partial derivatives of the 
optical flow are used to compute the motion of a camera [6], [9]. It has been claimed 
that it is sufficient to know the optical flow and both its first and second derivatives 
at a single point to uniquely determine the motion [9]. This is incorrect (except for 
a special case) [1]. Furthermore, noise in the measured optical flow is accentuated by 

differentiation. 

In the continuous approach, the whole optical flow field is used. A major shortcom¬ 
ing of both the local and differential approaches is that neither allows for errors in the 
optical flow data. This is why we choose the continuous approach and devise a least- 
squares technique to determine the motion of the camera from the measured optical 
flow. The proposed algorithm takes the abundance of available data into account and 
is robust enough to allow numerical implementation. 

2 . Technical Prerequisites 

In this section we review the equations describing the relation between the motion 
of a camera and the optical flow generated. We use essentially the same notation as 
[9]. A camera is assumed to move through a static environment. Let a coordinate 
system X,Y, Z be fixed with respect to the camera, with the Z-axis pointing along 
the optical axis. If we wish, we can think of the environment as moving in relation 
to this coordinate system. Any rigid body motion jean be resolved into two factors, 
a translation and a rotation. We will denote by T the translational component of 
the motion of the camera and by Co its angular velocity (see also Figure 1 which is 
redrawn from [9]). Let the instantaneous coordinates of a point P in the environment 

be (X,Y,Z). 
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Figure 1 . Coordinate Systems 


(Note that Z > 0 for points in front of the imaging system.) Let f he the vector 
(X,Y,Z) t , where T denotes the transpose of a vector, then the velocity of P with 
respect to the X,Y, Z coordinate system, is: 

V = — f — CSX?. ( 1 ) 

We define the components of T and Q as: 

f = (U,V,W) T cD = (A,B, C) T . (2) 

Thus we can rewrite (1) in component form: 

X' = — U — BZ + CY 

Y> = —y — CX + AZ (3) 

Z' = — W — ay + bx. 

where ' denotes differentiation with respect to time. 

The optical Sow at each point in the image plane is the instantaneous velocity of 
the brightness pattern at that point. Let ( x, y ) denote the coordinates of a point in the 
image plane (see Figure 1). Since we assume perspective projection between an object 
point P and the corresponding image point p, the coordinates of p are: 


mmsami 
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i. 





_X _Y 

x ~ z y ~z 

The optical flow, denoted by ( u,v ), at a point ( x,y ) is: 

u = x' v — y'. (5) 

Differentiating (4) with respect to time and using (3) we obtain the following equations 
for the optical flow: 



X' 

u ~ z 

XZ' _ 

£2 — V 

u 

z 

-B + Cy)- 

W 

- x[ - 

v Z 

- Ay + Bx ) 

(6) 

Y' 

v ~ z 

YZ' 

Z 2 ~ 

V 

z 

- Cx -(- A) - 

, w 

-y(~~ - 
Zj 

-Ay + Bx). 



We can write these equations in the form: 


u 


Ut “j - u r 


V 


v t 4 - v, 


p) 


where (u t ,v t ) denotes the translational component of the optical flow and (u r ,v r ) the 
rotational component: 


u t 


u. 


U + xW 


v t 


V + yW 


z z 

Axy — B(x 2 + 1 ) 4 Cy v r = A(y 2 -f 1) — Bxy — Cx. 


( 8 ) 


So far we have considered a single point P. To define the optical flow globally we 
assume that P lies on a surface defined by a function Z — Z[X,Y ) which is positive for 
all values of X and Y. With any surface and any motion of a camera we can therefore 
associate a certain optical flow and we say that the surface and the motion generate 
this optical flow. 

Optical flow, therefore, depends upon the six parameters of motion of the camera 
and upon the surface whose images are analyzed. Can all these unknowns be uniquely 
recaptured solely from optical flow? The answer is no. To see this, consider a surface 
S 2 which is a dilation by a factor k of a surface S\. Further, let two motions denoted by 
Mi and M 2 have the same rotational component and let their translational components 
be proportional to each other by the same factor k (we will say that Mi and M 2 are 
similar). Then the optical flow generated by Si and Mi is the same as the optical flow 
generated by S 2 and M 2 . This follows directly from the definition of optical flow (8). 
It is still an open question whether there are any other pairs of distinct surfaces and 
motions which generate the same optical flow. 

Determining the motion of a camera from optical flow is much easier if we are told 
that the motion is purely translational or purely rotational. In the next two sections 
we will deal with these two special cases. Then we shall analyze the case where no a 
priori assumptions about the motion of the camera are made. 





5 


3. Translational Case 


In this section we discuss the case where the motion of the camera is assumed to 
be purely translational. As before, let T — (U,V,W) be the velocity of the camera. 
Then the following equations hold (see (8)): 



—U + xW 
Z 



—V + yW 
Z 



3.1. Similar Surfaces and Similar Motions 


It will be shown next that if two flows generate the same optical flow, and we know 
that the motions are purely translational, then the two surfaces are similar and the two 
camera motions are similar. Let Z\ and Z 2 be two surfaces and let T'i — {Ux,Vx,Wx) T 
and ?2 = {U 2 , V 2 , W2) t define two different motions of a camera, such that Z \ and ?\ 
and Z 2 and T 2 generate the same optical flow: 


U\ -f- xWi —Vi -I - yW 1 

- v =-—-, 

Z x Z x 

U2 + xW 2 -V 2 + yW 2 

- v — - —- . 

Z2 %2 

Eliminating Zx, Z2, u and v from these equations we obtain: 

—Ux xWx — U2 xW% 

— Vi H - yWx — V2 -j- yWj 

We can rewrite this equation as: 

(—Ux + xWx )(—Vjj + J/k V2) = (— U2 + xW2)(—Vx + yWx), 


u — 

u = 


( 10 ) 

( 11 ) 



or: 

U\V 2 — xV 2 Wx — yUxW 2 + xyWxW 2 — U 2 Vx — xVxW 2 — yU 2 Wx + xyW 2 Wx- (14) 

Since we assumed that Zx and and Z 2 and T 2 generate the same optical flow, the 
above equation must hold for all x and y. Therefore the following equations have to 
hold: 

UxV 2 = U 2 Vx 

—V 2 Wx = —V 1 W 2 (15) 

-UxW 2 = -U 2 Wx. 


These equations can be rewritten as: 

Ux-Vx-Wx = U 2 :V 2 :W 2 






from which, it follows that Z 2 is a dilation of Z x . It is clear that the scaling factor 
between Z\ and Z 2 (or equivalently between and 2 2 ) cannot be recovered from the 
optical flow, regardless of the number of points at which the flow is known. By uniquely 

determining the motion of the camera, we will mean that the motion is uniquely 
determined up to a constant scaling factor. 


3.2. Least-Squares Formulation 

In general, the direction of the optical flow at two points in the image plane 
determine the motion of a camera in pure translation uniquely. There is a drawback 
however to utilizing so little of the available information. The optical flow we measure 
is corrupted by noise and it is desirable to develop a robust method which takes this 
into account. Thus we suggest using a least-squares method [4], [12] to determine the 
movement parameters and the surface (i.e., the best fit with respect to some norm). 

For the following we assume that the image plane is the rectangle xe[—w,w) and 
y[—h,h\. The same method applies if the image has some other shape. (In fact, it 
can be used on sub-images corresponding to individual objects in the case that the 
environment contains objects which may move relative to one another). Furthermore 
we have to assume that 1/Z is a bounded function and that the set of points where 1 jZ 
is discontinuous is of measure zero. This condition on l/Z assures us that all necessary 
integrations can be carried out. We wish to minimize the following expression: 


u 


W 


[(u 


+(. 


W 




(17) 


In this case then, we determine the best fit with respect to the ML 2 norm which is 
defined as: 

/ h nw 

/ [/(z>y)] 2 dxdy. (18) 

— hJ — w 

The steps in the least-squares method are as follows: First we determine that Z which 
minimizes the integrand of (17) at every point (x, y). Then we determine the values of 
U, V and W which minimize the integral (17). 

Let us introduce the following abbreviations: 


a = —U -f xW j3 = —V + yW. 


Note that the expected flow, given U, V and W is simply: 


a 



Then we can rewrite (17) as: 





( 20 ) 


h 


( 21 ) 
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Figure 2. Geometrical Interpretation 


We proceed now with the first step of our minimization method. Differentiating the, 
integrand of (17) with respect to Z and setting the resulting expression equal to zero 

yields: 


(« 


Therefore we can write Z as: 


^ \ ^ I / 


/M 


z’z 2 


o. 


( 22 ) 


a 2 -f- 0 2 
uni -I- vQ 


(23) 


This equation, by the way, imposes a constraint on U, V and W , since Z must be 
positive. We do not make use of this except to help us pick amongst two opposite 
solutions for the translational velocity later on. Note that now: 


a a uj3 — va ^ P _ u(3 — va 

U ~ Z ~ 15 a 2 + fP V Z~ a a 2 +p 2 

and we can therefore rewrite (17) as: 


(24) 



(ufi — va) 2 
a 2 ~T /? 2 


dxdy. 



It should be clear, by the way, that uniformly scaling U , V and IT does not change 
the value of (25). This is a reflection of the fact that we can determine the motion 
parameters only up to a constant factor. 

Before proceeding with the second step, we give a geometrical interpretation in 
Figure 2 of what we have so far. Suppose that the motion parameters U, V, and 





W are given. At any given point, say (xo>yo)> optical flow depends not only upon the 
motion parameters but also upon the value of Z at that point, Zq say. However, the 
direction of (u,v) does not depend upon Z 0 . The point (u,v) must lie along the line L 
in the uv-plane defined by the equation uf3 — vot = 0. Let the measured optical flow 
at ( x o,!/o) be denoted by (u m ,u m ), and let the closest point on the line L be {u a ,v a ). 
This corresponds to a particular Z a given by (23). The remaining error is the distance 
between the point (i£ m ,v m ) and the line L. The square of this distance is given by the 
integrand of (25). 

For the second step, we differentiate (25) with respect to U, V and W and set the 
resulting expressions equal to zero: 


h 


L 


>W 

w 
•w 


/ n pw 
~hJ —w 


(ya 


(3{u(3 — va)(ua -f- v/3) 
(a 2 + (3 2 ) 2 

a(u/3 — va)(uct -j- vf3) 
w (a 2 + /? 2 ) 2 

x/3)(ufi — uq:)(uq: + v/3) 


(a 2 -f- /3 2 ) 2 


dxdy — 0 


dxdy = 0 


dxdy = 0. 


Let us introduce the following abbreviation: 


K 


(u/? — va)(ua + v/3) 
( a 2 + y 0 2)2 


Then equations (26) can be rewritten as: 


(26) 


(27) 



l—V + y W)K)dxdy = 0 



-U + xW)K)dxdy = 0 



yU -f- xV)K]dxdy 




The sum of U times the first integral, V times the second integral, and W times the 
third integral is identically zero. Thus the three equations are linearly dependent. This 
is to be expected, for if: 


f{kU, kV, kW) — f(U, V, W), 
where / is a differentiable function and k a constant, then: 




0 . 


(30) 





g 


The result is also consistent with the fact that only two equations are needed, since the 
translational velocity can be determined only up to a constant factor. Unfortunately 
equations (28) are nonlinear in U, V and W and we are not able to show that they have 
a unique (up to a constant scaling factor) solution. 

3.3. Using a Different Norm 

There is a way, however, to devise a least-squares method which allows us to display 
a closed form solution for the motion parameters. Instead of minimizing (17), we will 
try to minimize the following expression: 

fh rw —£/ _j_ xW *2 I ( —V -f- yW . 2 1/ 2 I fl2\ J j (o-i\ 

/ / [(«-J ) +(u- z ) ](<* + P )dxdy (31) 

J — hJ—w " " 

obtained by multiplying the integrand of (17) by a 2 + /? 2 . Then we apply the same 
least-squares method as before to (31). When the measured optical flow is not corrupted 
by noise, both (31) and (17) can be made equal to zero by substituting the correct motion 
parameters. We thus obtain the same solution for the motion parameters whether we 
minimize (31) or (17). If the measured optical flow is not exact, then using expression 
(31) for our minimization, we obtain the best fit with respect not to the ML 2 norm, 
but to another norm which we call the ML a p norm: 


/ h rw 

/ [/(*, 2/)] V 2 + P 2 )dxdy. (32) 

-hJ —w 

What we have here is a minimization in which the error contributions are weighted, 
greater importance being given to points where the optical flow velocity is larger. This 
is most appropriate when the measurement of larger velocities is more accurate. 

Which norm gives the best results depends on the properties of the noise in the 
measured optical flow. The first norm is better suited to the sitation where the noise in 
the measurements is independent of the magnitude of the optical flow. Note also that 
if we really want the minimum with respect to the ML 2 norm, we can use the results 
of the minimization with respect to the ML a p norm as starting values in a numerical 
minimization. 

We discuss now our least-squares method in the case where the norm is chosen to 
be ML a p. First we determine Z by differentiating the integrand of (31) with respect 
to Z and setting the result equal to zero. We again get (22): 



a. a 
l)~Z 2 



l)£ 

z } z 2 



from which it follows that (23): 


Q 2 + ft 2 

ua -f- v/3 


( 34 ) 




So we want to minimize: 



(lift — va) 2 dxdy. 


Let us call this integral g(U, V, W), then, since: 



u0 — va = (vU 

— uV ) — [xv — yu)W, 

(36) 

we have: 



g{U, V, W) — aU 2 + bV 2 + cW 2 + 2 dUV + 2 eVW + 2 fWU, 

(37) 

where: 



ph pw 

a — 1 I v 2 dxdy 

J — hJ — w 

ph pw 

b— I / u 2 dxdy 

J — h<J —w 


ph pw 

c — I I (xv — yu) 2 dxdy 

J-hJ—w v ; 

ph pw 

d— — / / uvdxdy 

J — hJ —w 

(38) 

ph pw 

e — 1 u i xv — yu)dxdy 

J—hJ—w J 

ph pw 

f — — / / v{xv — yu)dxdy. 

J — h J —\d 



Now g(U, V, W) cannot be negative, and g{U, V, W) = 0 for U — V — W — 0. Thus 
a minimum can be found by inspection, but is not what we might have hoped for. In 
fact, to determine the translational velocity using our least-squares method we have to 
solve the following homogeneous equation for T: 

Gt = 0 (39) 

where G is the matrix: 

fa d f\ 

G — ( d b el. (40) 

\f e cJ 

Clearly (39) has a solution other then zero if and only if the determinant of G is zero. 
Then the three equations (39) are linearly dependent and f can be determined up to 
a constant factor. In general, however, as the data is corrupted by noise, g cannot be 
made equal to zero for non-zero translational velocity and so f — (0, 0, 0) T will be the 
only solution to (39). To see this in another way, note that g has the following form: 

g(kU, kV, kW) = k 2 g(U, V, W) (41) 

where k is a constant. Clearly g(U,V,W ) assumes its minimal value for U = V — 
W = 0 . 

What we are really interested in, is determining the direction of f which minimizes 
g, for a fixed length of f. Hence we impose the constraint that T be a unit vector. 
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If f is constrained to have unit magnitude, the minimum value of g is the smallest 
eigenvalue of the matrix G and the value of ? for which g assumes its minimum can be 
found by determining the eigenvector corresponding to this eigenvalue [8]. This follows 
from the observation that g is a quadratic form which can be written as: 


g{U, V, W) — T T GT. (42) 

Note that G is a positive semidefinite hermitian matrix as a > 0, b > 0, c > 0, ab > 
d 2 , he > e 2 and ca > f 2 . (The last three inequalities follow from the Cauchy-Schwarz 
inequality [7], [8]). Hence all eigenvalues are real and non-negative and are the solutions 
X of the third degree polynomial: 

X 3 — (a -f- b -f- c)X 2 -f- (ab -f- be -f- ca — d 2 — e 2 — / 2 )X -f- (43) 

( ae 2 -f- bf 2 + cd 2 — abc — 2 def) — 0. 

There is an explicit formula for the least positive root in terms of the real and imaginary 
parts of the roots of the quadratic resolvent of the cubic. In our case this gives us the 
desired smallest root, since the roots cannot be negative. For the sake of completeness, 
however, various pathological cases that might come up will be discussed next, even 

though they are of little practical interest. 

Note that X = 0 is an eigenvalue if and only if G is singular, that is, if the constant 
term in the polynomial (43) equals zero. In fact, if the determinant of G is zero one 
can find a velocity f which makes g zero. It follows from a theorem in calculus that 
this happens only when the optical flow is either not corrupted by noise at all or only 
at a few points. The theorem states that if the integral of the square of a bounded and 
continuous function is zero then the function itself is zero. Hence errors can only occur 
at points where the optical flow is discontinuous, and these are exactly the points where 
the surface defined by Z is discontinuous. (These are also the places where existing 
methods for computing the optical flow [5] are subject to large errors). 

It is impossible for exactly two eigenvalues to be zero, since this would imply that 
the coefficient of X in the polynomial (43) equalled zero, while the coefficient of X 2 did 
not. That in turn would imply that ab — d 2 , be = e 2 , and ca = / 2 , while a, b, and 
c are not all zero. For equality to hold in the Cauchy-Schwarz inequalities, however, 
u and v must both be proportional to xv — yu. This can only be true (for all x and 
y i n the image) if u=v— 0. But then all six integrals become zero and consequently 
all three eigenvalues are zero. This situation is of little interest, since it occurs only 
when the optical flow data is zero everywhere. Then the velocity is zero too. Once 
the smallest eigenvalue is known, it is straightforward to find the translational velocity 
which best matches the given data. To determine the eigenvector corresponding to an 
eigenvalue, say Xi, we have to solve the following set of linear equations: 





( 44 ) 


(a — \ t )U + dV + fW = 0 
dU + (b — \x)V + eW = 0 
fU + eV + {c — \ 1 )W = 0. 

As Xj is an eigenvalue, equations (44) are linearly dependent. Let us for a moment 
assume that all eigenvalues are distinct, that is, the rank of the matrix (G—\I), where 
I is the identity matrix, is two. Then we can use any pair of them to solve for U, V in 
terms of W say. There are three ways of doing this. For numerical accuracy we may 
add the three results to get the symmetrical forms: 



— /(* — Xi) - 

- d(c - 

“ Xi) -f- e(f -f- d * 

-e) 


v = (c — X x )(o — X x ) 

— d(c — Xi) - 

- e(a - 

~ ^i) “1" f(d + e 

-/) 

(45) 

W = (a — Xi)(6 — Xi) 

— e(a —Xi)- 

-}{b- 

~ Xi) -f~ d[e -f- f 

— d). 



Note that Xi will be very small, if the data is good, and one may wish to just 
approximate the exact solution by using the above equations with Xt set to zero. (Then 
there is no need to find the eigenvalue). In any case, the resulting velocity may now be 
normalized so that its^magnitude equals one. There is one remaining difficulty, arising 
from the fact that if T is a solution to our minimization problem, so is — T. Only one 
of these solution will correspond to positive values of Z in equation (34) however. This 
can be easily seen by evaluating (34) at some point in the image. The case where the 
two smallest eigenvalues are the same will be discussed in one of the next paragraphs. 

There is a simple geometrical interpretation of what we have done so far. To this 
end we consider the surface defined by g(U, V, W\ = ky here Jfc is a constant. Note 
that we can always find a new coordinate system U, V, W in which g(U, V, W) can be 
written as: 

U -j- \%V -T X 3 W — k (46) 

where X,' for i = 1, 2, 3 are the three eigenvalues of the quadratic form. If the eigenvalues 
are all non-zero, the surface g(U, V, W) = k is an ellipsoid with three orthogonal semi¬ 
axes of length \/kJ\i. We are particularly interested in the case where the constant 
k is the smallest eigenvalue. Then all three semi-axes have lengths less than or equal 
to one. Hence the ellipsoid lies within the unit sphere. If the two smallest eigenvalues 
are distinct, the unit sphere touches the ellipsoid in two places, corresponding to the 
largest axis. If the two smaller eigenvalues happen to be the same, however, the unit 
sphere touches the ellipsoid along a circle and as a result all the velocity vectors lying 
in a plane spanned by two eigenvectors give equally low errors. Finally, if all three 
eigenvalues are equal, no direction for T is preferred, since the ellipsoid becomes the 
unit sphere. 

The case where exactly one eigenvalue is zero also has a simple geometrical inter¬ 
pretation. The surface defined by g(U, V, W) = 0 is a straight line, which can be seen 
easily from the following equation: 




13 


written for the case when X 3 is zero. (Remember that Xi and X 2 are both positive.) 
Clearly the unit sphere intersects this line in exactly two points, one of them cor¬ 
responding to positive values for Z in equation (34). 

The method which we just described can be easily implemented. To this end, the 
problem can be discretized. An expression similar to (31) can be derived where the 
integrals are approximated by sums. Our minimization method can then be applied 
to these sums. The resulting equations are similar to ones described in this section, 
with summation replacing integration. We implemented the resulting algorithm and 
tested it using synthetic data including additive noise. Thes results agreed with our 
expectations. 

One can use the ratio of the biggest to the smallest eigenvalue, the so called 
condition number [13], as a measure of confidence in the computed velocity. The result 
is very sensitive to errors in the measurements unless this ratio is much bigger than 
one. 

Curiously, the same error integral as (35) is obtained in the case where the MLz U v 
norm is used: 

/ h pw 

/ [/(x, y)Z(x, y)] 2 [u 2 + v 2 )dxdy. (48) 

-hJ —w 

We can arrive at a similar solution by multiplying the integrand in (17) by Z 2 instead 
of by a 2 + ft 2 - In that case the minimization is carried out with respect to the MLz 
norm defined by: 

/ 1 h pw 

/ [f{x,y)Z{x,y)] 2 dxdy. (49) 

-hJ —w 

Here optical flow velocities for points which are further away are weighted more heavily. 
This is most appropriate when the measurement of larger velocities is less accurate. 
We end up with a quadratic form similar to g, only the integrals for the six constants 
corresponding to a, b, c, d, e, and / are a bit more complicated. Curiously they only 
depend on the direction of the optical flow at each point, not its magnitude. 

Also, other constraints could be used. If we insist on U 2 V 2 — 1, for example, we 
obtain a quadratic instead of a cubic equation, and if we use W = 1 , a linear equation 
only need to be solved. The disadvantage of these approaches is that the result is 
sensitive to the orientation of the coordinate axes. Clearly, in the case of exact data, 
we get the right solution using any of the constraints mentioned above. 


3.4. Using a Different Constraint 

The minimization scheme discussed in the previous section gives us a unique 
solution in most cases for the velocity vector T. Here we propose a slightly dtffere 
approach which always gives us a unique solution. Note that applying the first step 
our nunirnizafior! method gives us a constraint between the values of Z, the veloci 




vector and the optical flow at every point. We can in addition assume that Z — Z 0 is 

known at a particular point, say (xo>yo)- Using the MLz U v norm in our scheme, we 
want to minimize: 



uZ — (—U + xW)} 2 + 


[vZ - {—V 4 #)] 2 (u 2 + v 2 )dxdy. 



Differentiating (50) with respect to Z, and setting the resulting expression equal to 
we obtain: 

— uoc —J- v/3 

Zj ' ==z —--—r-. 

U 2 + V 2 


zero, 

(51) 


Thus we propose to solve the minimization scheme under the following constraint: 


Zo(uo 4 v o ) — (“o a 4 v oP) — 0 (52) 

where uq and uq denote the components of the optical flow measured at (xoil/o)- The 
error integral (50) becomes after substituting (51): 



(u/? — va) 2 dxdy 



which is the same as (35) and is denoted by g{U, V, W) (37). Thus we want to minimize: 


g{U, V , W ) 4- 2X[Zq(uq 4" u o) — (wo<* 4- ^o/?)] 



where X denotes a Lagrangian multiplier. To determine U, V and W the following 
linear equations obtained by differentiating (54) with respect to U, V, W and X have to 
be solved: 


aU 4 dV 4 fW 4 \u 0 = 0 

dU 4 bV 4- eW 4 Xu 0 = 0 (55) 

fU 4 eU 4- cW - X(z 0 u 0 4 Vovo) = 0 

Wo U 4 v 0 V — (x 0 uq 4- y 0 v 0 )W = — Z 0 {u 2 0 4- v 2 ). 

These equations can be written in the form: 


G x ?x = P (56) 

where T\ — (U,V,W,\) T and P — (0,0,0,—^o( w o "U v o)) T - Let the determinant of 
G\ be Ao: 

A 0 ={d 2 — ab)(x 0 u 0 4 y 0 v 0 ) 2 4- (e 2 — 6c)uj§ 4- (/ 2 — ac)v 2 0 4 (57) 

2\{de — bf)u 0 (x 0 u 0 4 2/o v o) 4 {df — (ie)v 0 (x 0 u 0 4 Vo^o) 4- (cd — ef)u 0 v 0 ]. 
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Assuming that Ad ^ 0 m can easily determine T\ from (55): 




(58) 


Introducing the following abbreviation: 


Zq( u q 4~ «q) 
Ao 


(59) 


we can give these formulae for 7V: 


K[uo(bc — e 2 ) -f- vq (ef — cd) -f- (zo^o -j- yoVo)(bf — ^ e )] 
K\uo{ef — cd) -f- vq (ac — f 2 ) -f- (xo u o yo u o)( ae — df)\ 
K[uo(de — bf) -f- vq (df — ae) -f- (^o^o + 2/o u o)(^ 2 a ^)} 

K[ae 2 + d 2 + bf 2 — abc — 2 def]. 


(60) 


The disadvantage of this approach is that the result depends upon the values of the 
optical flow at a single point. To circumvent this problem we propose to determine 
average values for U , V and W in the following manner. First note that we are only 
interested in the ratios of UjW and VfW which obviously do not depend upon the 
(unknown) value for Zq. Equivalently we could determine the value for K from the 
condition that f should have unit length. Hence we can determine values for U,V 
and W which depend only upon the values of the optical flow at a single point and the 
coefficient in the matrix G. If we want to remove the dependence of the result on the 
data at a single point, we can simply average the values obtained for all image points. 

In the case where the data is exact, i.e., where the determinant of G, denoted 
by detG, vanishes, the solution for T is the same one as obtained using no constraint 
in our minimizations scheme. To see this just observe that in that case X = 0 as 
X = —KdetG). In the case where A 0 = 0, equations (55) have a solution only when 
detG = 0. We do not have to be concerned with the case where A 0 = 0 but detG yZ 0 
as we can argue that equations (55) always have to have a solution. Note that our 
method is based on the condition that Z is a certain function (51) of U,V ,W. Hence 
(52) cannot impose a constraint which would be impossible to satisfy. 

4. Rotational Case 

Suppose now that the motion of the camera is purely rotational. In order to 
determine the motion from optical flow we again use a least-squares algorithm with the 
ML 2 norm described in the previous section. Recall that in this case the optical flow 
is (see (8)): 


Axy — B(x 2 + 1) + Gy 
A{y 2 + 1) — Bxy - Cx. 


( 61 ) 
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We will show now in an analogous fashion to section 3.1 that two different rotations, 
say S>i = (Ai,Bi,Ci) t and u >2 = (A 2 , B 2 , C 2 ) T , cannot generate the same optical flow. 
Assuming the converse, the following equations have to hold for all values of x and y: 


A\xy -j~ 1 } -f- C\y — A2xy B2(x 2 -j- 1 ) -j - C*2 y 

A i[y H - l) B\xy C\x — A2(y 2 1 ) — B2xy — C2X 


(62) 


from which we can immediately deduce that uii — cD 2 - 

In general, the direction of the optical flow at two points and its magnitude at one 
point determine the motion of a camera in pure rotation uniquely. We choose instead 
to minimize the following expression: 


ph pw 

/ J K u ~ Wr ) 2 + ( v — Vr ) 2 ] 

J —hJ —w 


dxdy. 


(63) 


As the motion is purely rotational, the optical flow does not depend upon the distance 
to the surface and therefore we may omit the first step in our method. Thus we im¬ 
mediately differentiate (63) with respect to A, B and C and set the resulting expressions 
equal to zero: 


rh pw 

/ / [( w — Ur )xy -f (w — 1>, 

J — hJ —w 
ph pw 

/ ,/ [(u — U r )(x 2 + 1) -f (1 

J - hj —~W 


u r )xy -f (v — v r )(y 2 -f 1 )]dxdy 


ph pw 

[(u — u r )y — {v 

J —nJ — W 


hJ —w 


v r )xy]dxdy 
- v r )x]dxdy 


(64) 


Rewriting the above equations we obtain: 


nh nw ph pyj 

J_ h J_ w l ux y + v (y 2 + 1 ))dxdy = J J \u r xy + v r (y 2 -f 1 )]dxdy 


h nw 


f h r 

/ / [u(a; 2 -j- 1) + vxy]dxdy 

J — hJ —w 


hJ —w 


/ h pw 

J-„ luy 


p h pw 

I / [u r (z 2 + 1) + v r xy]dxdy 

J — h,J —W 
nh nw 


(65) 


hJ —w 


vx]dxdy 


/ n pw 
-hi —w 


u r y — v r x]dxdy 


and expanding these equations yields: 


aA -f dB + fC —k 
dA-\-bB + eC =7 
fA -f- eB -f- cC =fh, 


( 66 ) 
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where: 


ph pw 

/ / [x 2 y 2 + (y 2 + l) 2 ]dxdy 

J — hJ—w 


ph pw 

/ / [(z 2 + i) 2 + x 2 V 2 ]dxdy 

J — hJ —W 


-hJ —w 
h pw 


pn pw 

/ / [x 2 + y 2 ]dxdy 

J —frj - W 


(67) 


hJ —w 
ph pw 


pn pw 

/ / [xy(x 2 + V 2 + 2 ))dxdy 

J — hJ —w 


/ h pw 

/ ydxdy 
-hJ —w 


/ h pw 

I xdxdy, 
-hJ —w 


and: 


ph pw 

/ / [uxy + v(y 2 -f l)]dxdy 

J - fiJ - yj 


ph pw 

/ / [u(x 2 + 1) -f vxy]dxdy 

J - }%J - yj 


( 68 ) 


ph pw 

LJ-J uy 


vx]dxdy. 


If we call the coefficient matrix in (66) M and the column vector on the right-hand side 
71, then we have: 

MQ = 71. ( 69 ) 

Thus, provided the matrix M is non-singular, we can compute the rotation as follows: 


M" 1 n. 


(70) 


It is easy to see that the matrix M is non-singular in the special case of a rectangular 
image plane since then: 


4wh{ ^ + H|! + 1} + 0 ™. 


( 71 ) 




So in the case of a rectangular image plane, the matrix is diagonal, which makes it 
particularly easy to compute its inverse. In fact, the matrix is diagonal if the image 
plane is symmetrical about the z-axis and the y-axis. As the extent of the image plane 
decreases, however,^the matrix M becomes ill conditioned. That is inaccuracies in the 
three integrals (fc, /, and ttt.) computed from the observed flow are greatly magnified. 
This makes sense since we cannot expect to accurately determine the component of 
rotation about the optical axis when observations are confined to a small cone of 
directions about the optical axis. Again, in our numerical implementation of the 
algorithm the integrals in (67) can be approximated by sums. 


5. General Motion 

We would like now to apply a least-squares algorithm to determine the motion of a 
camera from optical flow given no a priori assumptions about the motion. It is plain that 
a least-squares method is easiest to use when the resulting equations are linear in all the 
motion parameters. Unfortunately, there exists no norm which will allow us to achieve 
this goal. We did find a norm, however, which resulted in equations that are linear in 
some of the unknowns and quadratic in the others. We again attack the minimization 
problem using the ML q/3 norm under the constraint that U 2 -f V 2 + W 2 = 1. The 
ensueing equations are polynomials in the unknowns U, V,W,A,B and C and can be 
solved by a standard iteration method like Newton’s method or Bairstow’s method [12] 

or by an interpolation scheme like regula falsi [12]. The expression we wish to minimize 
is: 

/ h pw 

_ h J_ w ^ u ~ + u r )] 2 + W — {-g -f- Vr)} 2 }(* 2 + /3 2 )dxdy. (72) 

The first step is to differentiate the integrand of (72) with respect to Z and set the 
resulting expression equal to zero: 



q 2 -f (3 2 

(u — u r )ot -f {v — v r )/3' 



We introduce the Langrangian multiplier X as before and attempt to minimize: 



u r )/3 — (v — 


v r )ct) 2 dxdy -f- X(I7 2 -f- V 2 -f- W 2 




The equations we have to solve to determine the motion parameters are obtained by 
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differentiation: 


(75) 


Note that the first three of these equations are linear in A, B and C from which these 
parameters can be determined uniquely in terms of U, V and W. Then we can determine 
U,V and W from the last four equations by a numerical method. To this end, the 
problem can be discretized and equations analogous to (75) derived, where summation 
of the appropriate expressions is used instead of integration. 

6. Summary 

Our objective was to devise a method for determining the motion of a camera from 
optical flow which allows for noise in the measured data. The least-squares method 
which we proposed in this paper meets our goal and is also suitable for numerical 
implementation. An important application of our results is in passive navigation. Here 
the path and instantaneous altitude of a vehicle is to be determined from information 
gleaned about the environment without the emission of sampling radiation from the 
vehicle. 
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